function m = cmoment(bw, x, y, si)

% convert bw image to indices
[i j] = find(bw == 1);

% subtract mean
iMS = i - mean(i);
jMS = j - mean(j);

% power the mean-subtracted components to the appropriate order
iMSP = iMS .^ y;
jMSP = jMS .^ x;

% compute moment
m = sum(iMSP .* jMSP);

if si
   m = m / (cmoment(bw, 0, 0, 0) ^ (((x + y) / 2) + 1)); % normalize for scale
end

end